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ABSTRACT 



The multi-frequency capability of the Planck satellite provides information both on the integrated history of star formation (via the cosmic infrared 
background, or CIB) and on the distribution of dark matter (via the lensing effect on the cosmic microwave background, or CMB). The conjunction 
of these two unique probes allows us to measure directly the connection between dark and luminous matter in the high redshift (1 < z < 3) 
Universe. We use a three-point statistic optimized to detect the correlation between these two tracers. Following a thorough discussion of possible 
ff) contaminants and a suite of consistency tests, using lens reconstructions at 100, 143 and 217 GHz and CIB measurements at 100-857 GHz, we 
report the first detection of the correlation between the CIB and CMB lensing. The well matched redshift distribution of these two signals leads 
£f~) to a detection significance with a peak value of 42 a at 545 GHz and a correlation as high as 80 % across these two tracers. Our full set of multi- 
t— H frequency measurements (both CIB auto- and CIB-lensing cross-spectra) are consistent with a simple halo-based model, with a characteristic mass 
scale for the halos hosting CIB sources of log 10 (M/M Q ) = 10.5 ± 0.6. Leveraging the frequency dependence of our signal, we isolate the high 
redshift contribution to the CIB, and constrain the star formation rate (SFR) density at z > 1. We measure directly the SFR density with around 
2 cr significance for three redshift bins between z = 1 and 7, thus opening a new window into the study of the formation of stars at early times. 
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1. Introduction 

This paper, one of a set associated with the 2013 release of 
data from the Planck 1 mission (Planck Collaboration I 2013), 



* Corresponding author: Olivier Dore <olivier .p . dore@jpl . 
nasa.gov> 

1 Planck (http://www.esa.int/Planck) is a project of the 
European Space Agency (ESA) with instruments provided by two sci- 
entific consortia funded by ESA member states (in particular the lead 



presents a first detection of a strong correlation between the in- 
frared background anisotropies and a lensing-derived projected 
mass map. The broad frequency coverage of the Planck satel- 
lite provides two important probes of the high redshift Universe. 
In the central frequency bands of Planck (70, 100, 143, and 
217 GHz), cosmic microwave background (CMB) fluctuations 



countries France and Italy), with contributions from NASA (USA) and 
telescope reflectors provided by a collaboration between ESA and a sci- 
entific consortium led and funded by Denmark. 
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dominate over most of the sky. Gravitational lensing by large- 
scale structure produces small shear and magnification effects in 
the observed fluctuations, which can be exploited to reconstruct 
an integrated measure of the gravitational potential along the line 
of sight Okamoto & Hu (2003). This "CMB lensing potential" 
is sourced primarily by dark matter halos located at 1 < z < 3, 
halfway between ourselves and the last scattering surface (see 
Blandford & Jaroszynski 1981; Blanchard & Schneider 1987, or 
Lewis & Challinor 2006 for a review). In the upper frequency 
bands (353, 545, and 857 GHz), the dominant extragalactic sig- 
nal is not the CMB, but the cosmic infrared background (CIB), 
composed of redshifted thermal radiation from UV-heated dust, 
enshrouding young stars. The CIB contains much of the energy 
from processes involved in structure formation. According to 
current models, the dusty star-forming galaxies (DSFGs), which 
form the CIB have a redshift distribution peaked between z ~ 1 
and z ~ 2, and tend to live in 10 n -10 13 M Q dark matter halos 
(see, e.g., Bethermin et al. 2012, and references therein). 

As first pointed out by Song et al. (2003), the halo mass and 
redshift dependence of the CMB lensing potential and the CIB 
fluctuations are well matched, and as such a significant correla- 
tion between the two is expected. This point is illustrated quan- 
titatively in Fig. 1, where we plot estimates for the redshift- and 
mass- kernels of the two tracers. In this paper we report on the 
first detection of this correlation. 

Measurements of both CMB lensing and CIB fluctuations 
are currently undergoing a period of rapid development. While 
the CIB mean was first detected using the FIRAS and DIRBE 
instruments aboard COBE (Puget et al. 1996; Fixsen et al. 1998; 
Hauser et al. 1998), CIB fluctuations were later detected by 
the Spitzer Space Telescope (Lagache et al. 2007) and by the 
BLAST balloon experiment (Viero et al. 2009) and the Herschel 
Space Observatory (Amblard et al. 2011; Viero et al. 2012), 
as well as the new generation of CMB experiments, includ- 
ing Planck, which have extended these measurements to longer 
wavelengths (Hall et al. 2010; Dunkley et al. 2011; Planck 
Collaboration XVIII 2011; Reichardt et al. 2012). The Planck 
early results paper: Planck Collaboration XVIII (2011) (hence- 
forth referred to as PER) presented measurements of the angu- 
lar power spectra of CIB anisotropies from arc-minute to degree 
scales at 217, 353, 545, and 857 GHz, establishing Planck as a 
potent probe of the clustering of the CIB, both in the linear and 
non-linear regimes. A substantial extension of PER is presented 
in a companion paper to this work (Planck Collaboration 2013, 
henceforth referred to as PIR). 

The CMB lensing potential, on the other hand, which was 
first detected statistically through cross-correlation with galaxy 
surveys (Smith et al. 2007; Hirata et al. 2008, and more recently 
Bleem et al. 2012; Sherwin et al. 2012), has now been observed 
directly in CMB maps by the Atacama Cosmology Telescope 
and the South Pole Telescope (Das et al. 201 1; van Engelen et al. 
2012). 

Planck's, frequency coverage, sensitivity and survey area, al- 
low high signal-to-noise measurements of both the CIB and the 
CMB lensing potential. Accompanying the release of this pa- 
per, Planck Collaboration XVII (2013) reports the first measure- 
ment and characterisation of the CMB lensing potential with the 
Planck data, which has several times more statistical power than 
previous measurements, over a large fraction (approximately 
70% of the sky). We will use this measurement of the lensing 
potential in cross-correlation with measurements of the CIB in 
the PlanckHFl bands to make the first detection of the lensing- 
infrared background correlation. In addition to our measure- 
ment, we discuss the implications for models of the CIB fluc- 
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Fig. 1. Redshift- and mass- integrand for the CIB and CMB lens- 
ing potential power spectra at ( — 500, calculated using the 
CIB halo model of Planck Collaboration XVIII (2011), evalu- 
ated at 217 GHz. The good match between the redshift and halo 
mass distributions leads to an expected correlation up to 80%. 
The sharper features in the CIB kernel are artefacts from the 
Bethermin et al. (2012) model. We note that the low mass, high 
z behavior of our calculation is limited by the accuracy of the 
mass function we use (Tinker & Wetzel 2010). All of our mass 
integrals use M,„,„ = 10 5 M Q . 

tuations. The outline of this paper is as follows. In Sect. 2 we 
describe the data we will use, followed by a description of our 
pipeline for correlating the CIB and lensing signals in Sect. 3. 
Our main result is presented in Sect. 4, with a description of our 
error budget, consistency tests and an array of systematic tests in 
Sect. 5. We discuss the implications of the measured correlation 
for CIB modelling in Sect. 6. 

2. Data sets 

2. 1 . Planck maps 

Planck (Tauber et al. 2010; Planck Collaboration I 2011) is the 
third generation space mission to measure the anisotropy of the 
CMB. It observes the sky with high sensitivity in nine frequency 
bands covering 30-857 GHz at an angular resolution from 31' to 
5'. The Low Frequency Instrument (LFI; Mandolesi et al. 2010; 
Bersanelli et al. 2010; Mennella et al. 2011) covers the 30, 44, 
and 70 GHz bands with radiometers that incorporate amplifiers 
cooled to 20 K. The High Frequency Instrument (HFI; Lamarre 
et al. 2010; Planck HFI Core Team 2011a) covers the 100, 143, 
217, 353, 545, and 857 GHz bands with bolometers cooled to 
0. 1 K. Polarization is measured in all but the highest two bands 
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(Leahy et al. 2010; Rosset et al. 2010). A combination of radia- 
tive cooling and three mechanical coolers produces the temper- 
atures needed for the detectors and optics (Planck Collaboration 
II 2011). Two data processing centres (DPCs) check and cali- 
brate the data and make maps of the sky (Planck HFI Core Team 
2011b; Zacchei et al. 2011). Planck's sensitivity, angular reso- 
lution, and frequency coverage make it a powerful instrument 
for Galactic and extragalactic astrophysics as well as cosmol- 
ogy. Early astrophysics results are given in Planck Collaboration 
VIII-XXVI 201 1, based on data taken between 13 August 2009 
and 7 June 2010. Intermediate astrophysics results are now be- 
ing presented in a series of papers based on data taken between 
13 August 2009 and 27 November 2010. This paper uses data 
corresponding to the second Planck data release, with data ac- 
quired in the period up to 27 November 2010 and undergoing 
improved processing. 

We use the Planck HFI temperature maps at all six frequen- 
cies, i.e., 100, 143, 217, 353, 545, and 857 GHz. The maps at 
each frequency were created using almost three full-sky sur- 
veys. Here we give an overview of the HFI map-making pro- 
cess with additional details given in Planck HFI Core Team 
(201 lb); Planck Collaboration VI (2013). The data are organized 
as time-ordered information, hereafter TOI. The attitude of the 
satellite as a function of time is provided by two star trackers 
on the spacecraft. The pointing for each bolometer is computed 
by combining the attitude with the location of the bolometer in 
the focal plane, as determined by planet observations. The raw 
bolometer TOI for each channel is first processed to produce 
cleaned timelines and to set flags that mark bad data (for ex- 
ample data immediately following a cosmic ray strike on the de- 
tector). This TOI processing includes: (1) signal demodulation 
and filtering; (2) deglitching, which flags the strong part of any 
glitch and subtracts the tails; (3) conversion from instrumental 
units (volts) to physical units (watts of absorbed power, after a 
correction for the weak non-linearity of the response); (4) de- 
correlation of thermal stage fluctuations; (5) removal of the sys- 
tematic effects induced by 4 K cooler mechanical vibrations; and 
(6) deconvolution of the bolometer time response. Focal plane 
reconstruction and beam shape estimation is made using obser- 
vations of Mars. The simplest description of the beams, an el- 
liptical Gaussian, leads to full-width at half-maximum (FWHM) 
values of 9.65, 7.25, 4.99, 4.82, 4.68 and 4.33 'as given in Table 
4 of Planck Collaboration VI (2013). As explained in this paper, 
the inter-calibration accuracy between channels is better than the 
absolute calibration. This leads us to adopt conservative abso- 
lute calibration uncertainties of 0.64, 0.53, 0.69, 2.53, 10., 10. % 
at 100, 143 217, 353, 545 and 857 GHz respectively. We con- 
vert between emissivities given in MJy sr _1 (with the photomet- 
ric convention vl v = constant) and temperatures in juK, using 
the measured bandpass filters (see PER and PIR for details). 

For the sake of consistency testing (presented in particular in 
Sect. 5), we will sometimes use temperature maps where only a 
fraction of the TOI is used to generate the sky map. In particular, 
throughout this paper we use the terminology "half-ring" (HR) 
maps to refer to maps made using the first and second half of the 
stable pointing period, "survey" for individual full-sky survey 
maps (note that the third survey is incomplete for the particu- 
lar data release used in the intermediate papers), and "detset" 
for maps made using two independent sets of detectors per fre- 
quency (for details see Planck HFI Core Team 201 lb). 

We create three masks to exclude regions with bright 
Galactic and extragalactic foreground emission. The first mask 
accounts for diffuse Galactic emission as observed in the Planck 
data. To allow us to test for the effects of residual Galactic 




Fig. 2. Combined Galactic, point-source and H i mask with sky 
fractions 16, 30 and 43 %. 



emission on our results we create three different versions of 
this mask, each with a different masked area, such that 20, 40 
or 60 % of the sky is unmasked. Each version of this mask 
is created directly from the Planck 353 GHz map, from which 
we remove the CMB using the 143 GHz channel as a CMB 
template before smoothing by a Gaussian with FWHM of 5°. 
The map is then thresholded such that the mask has the re- 
quired sky fraction. Although the Galactic emission is stronger 
at 857 GHz, we expect the 353 GHz mask to better trace dust 
emission at the lower frequencies we use. The mask therefore ac- 
counts for Galactic dust and Galactic CO emission as explained 
in Planck Collaboration XII (2013). We will not worry about 
synchrotron emission, which is important at low frequencies, 
since its contribution at 100 GHz and at high Galactic latitudes 
is small, and, as with the dust component, will be uncorrelated 
with the lensing potential. The second mask covers bright point 
sources. This mask is created using algorithms tailored to de- 
tect point sources in the Planck data and is optimized for each 
frequency, as detailed in Planck Collaboration VII (2011) and 
Planck Collaboration (2011). The third mask is designed to re- 
move extended high-latitude Galactic dust emission ("cirrus"), 
as traced by external Hi data, as we will describe in Sect. 2.2.1. 
While the first two masks are described in Planck Collaboration 
XII (2013), the latter is specific to our cross-correlation analy- 
sis, as it provides a method to reduce the large-scale noise in our 
measurement, and the 3-point nature of our estimate ensures that 
it will not introduce a bias (although we test for this in Sect. 5). 
Ultimately, when we combine the three masks we obtain an ef- 
fective sky fraction of 16, 30 and 43 % for the 20, 40 and 60 % 
Galactic masks, respectively. 

2.2. External data sets 
2.2.1. H i maps 

We use measurements of 21 -cm emission from Galactic neutral 
hydrogen (H i) as a cirrus monitor. Outside of our Galactic and 
point source masks we use the H i data to construct a template 
of the dust emission in regions where the Hi column density 
is low (less than N^i < 2 x 10 20 cm" 2 ), and we mask regions 
where it is high, since in these regions the Hi and dust emis- 
sion are not well correlated (Boulanger et al. 1996; Boulanger & 
Perault 1988, PER). The masking procedure that we use is de- 
scribed in detail in Planck Collaboration XXIV (2011). It con- 
sists of subtracting the H i dust template from the Planck tem- 
perature map at 857 GHz and calculating the skewness of the 
residuals in 5 deg 2 regions. If the skewness is larger than a given 
value then the region is masked. This is an improvement over 
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the usual cut-off in H i column density. We use the latest release 
from the Leiden/ Argentina/Bonn (LAB) survey (Kalberla et al. 
2005), which consists of the Leiden/Dwingeloo Survey (LDS) 
(Hartmann & Burton 1997) north of -30° declination, com- 
bined with the Instituto Argentino de Radioastronomia Survey 
(Arnal et al. 2000; Bajaja et al. 2005) south of -25° declina- 
tion. The angular resolution of the combined map is approx- 
imately 0.6° FWHM. The LAB Survey is the most sensitive 
Milky Way H i survey to date, with the greatest coverage both 
spatially and kinematic ally. We make use of projections of the 
LAB maps onto =512 HEALPix maps performed by Land 
& Slosar (2007) and made available through the LAMBDA web- 
site 1 . The local standard of rest velocity coverage spans the in- 
terval -450kms~' to +400kms~', at a resolution of 1.3 kms -1 , 
with an rms brightness-temperature noise of 0.07-0.09 K, and 
with additional errors due to defects in the correction for stray 
radiation that are less than 20^1-0 mK for most of the data. 

2.2.2. IRIS//RAS maps 

As a consistency test we will use an additional tracer of the CIB 
that derives from re-processed IRAS maps at 60 and 100 fim. 
This new generation of IRAS maps, known as IRIS, benefits from 
improved zodiacal light subtraction, a calibration and zero level 
compatible with DIRBE, and an improved de-striping procedure 
(Miville-Deschenes & Lagache 2005). IRAS made two full-sky 
maps (HCON-1 and HCON-2), as well as a final map that cov- 
ers 75 % of the sky (HCON-3). The three maps had identical 
processing that included deglitching, checking of the zero-level 
stability, visual examination for glitches and artifacts, and zodia- 
cal light removal. The three HCONs were then co-added, taking 
into account the inhomogeneous sky coverage maps, to gener- 
ate the average map (HCON-0). Note that the Finkbeiner et al. 
(1999) maps are also constructed from the IRAS 100 //m data, 
and as such we will not investigate their cross-correlation prop- 
erties since the IRIS map contains the same information. For 
simplicity we will assume that the effective IRIS beam is uni- 
form across the sky and described by a Gaussian with FHWM of 
4.3'. 

3. Cross-correlation formalism and implementation 

We now describe our statistical formalism and its implementa- 
tion, with additional technical details given in the appendices. 
Our analysis consists of cross-correlating a full-sky reconstruc- 
tion of the CMB lensing potential with a temperature map. 

3. 1 . Reconstructing the CMB lensing potential 

The CMB is lensed by the gravitational potential of all matter 
along the photon trajectory from the last scattering surface to us. 
The lensed CMB is a remapping of the unlensed CMB with the 
lensed temperature equal to 0(« ) = ©(« + V0), where &(n ) 
is the unlensed CMB temperature and (p is the lensing poten- 
tial. We use the methodology described in Planck Collaboration 
XVII (2013) to obtain estimates (piM of the lensing potential 
in harmonic space, using the standard Okamoto & Hu (2003) 
quadratic estimator. 

Complete details on the lens reconstruction procedure, which 
we use are given in Planck Collaboration XVII (2013), although 
we review it briefly in point form here. Our estimates of (f> are 
obtained by the following set of steps: 

1 http : //lambda . gsf c . nasa . gov/product/f oreground/ 



1 . Inverse variance filter the CMB map. 

2. Use the filtered CMB map as the input to a quadratic lensing 
estimator, which is designed to extract the off-diagonal con- 
tributions to the CMB covariance matrix induced by lensing. 

3. Subtract a "mean-field bias", which corrects for known non- 
lensing contributions to the covariance matrix, including in- 
strumental noise inhomogeneity, beam asymmetry, and the 
Galaxy+point source mask. 

The output from this pipeline is an estimate of the lensing po- 
tential in harmonic space (/>lm and an associated noise power 
spectrum 2V£ , which we use to weight our cross-correlation es- 
timates. We also produce a set of simulated lens reconstruction, 
which we use to establish our statistical error bars. 

Our nominal lens reconstructions use the 143 GHz channel, 
however there is almost equivalent power to measure lensing us- 
ing the 217 GHz channel. Combining both channels would re- 
duce the noise power spectrum of our lens reconstruction by ap- 
proximately 25 %, compared with using either individually (the 
improvement is significantly less than 50 % because a signifi- 
cant portion of the lens reconstruction noise is due to the finite 
number of CMB modes, which we are able to observe, and is 
correlated between the two channels). We choose to focus on 
143 GHz here because it is significantly less susceptible to CIB 
contamination. We will use lens reconstructions based on the 
100 and 217 GHz data for consistency tests. 

3.2. Decreasing the foreground noise 

An important source of noise (but, as we will explain below, 
not bias) in our cross-correlation measurement is Galactic fore- 
ground emission. Dust emission is the dominant Galactic com- 
ponent at HFI frequencies above 217 GHz (see Sect. 5.1 for a 
quantitative discussion). In order to reduce the Galactic dust 
emission we create a dust template and subtract it from the tem- 
perature maps described in Sect. 2.1. At 100 and 143 GHz the 
CMB signal is significantly brighter than the dust emission out- 
side the Galactic mask. We therefore do not create and subtract 
a dust template at these frequencies. Note that while we could 
use other frequency maps to trace the CMB and remove it, to 
quantify the non-negligible amount of CIB that would be re- 
moved this way is not easy given the uncertainties in the cross- 
frequency CIB correlation structure. 

We rely on the well documented (but complex) correlation 
between Galactic Hi and dust (e.g., Boulanger & Perault 1988; 
Boulanger et al. 1996; Lagache et al. 1998, PER) to reduce the 
contamination by subtracting the H i-correlated dust component. 
As was performed in PER, we split the Hi map into two ve- 
locity components: a low -velocity local component (LC) typical 
of high-latitude H i emission, and a component of intermediate- 
velocity clouds (IVC). We found that the inclusion of a high- 
velocity component makes a negligible difference to the dust- 
cleaned map. Unlike the dedicated high-resolution H i observa- 
tions used in PER and PIR that only partially cover the sky, here 
we use the full-sky, low resolution LAB survey introduced in 
Sect. 2.2.1 as our Hi tracer. Although the resolution of this sur- 
vey is lower than the Planck resolution, it allows us to perform 
dust cleaning on large scales, where our cross-correlation mea- 
surement has high signal-to-noise ratio. The emissivity of the 
dust varies across the sky, and so the correlation between the 
dust and Hi emission is expected to vary. To account for this 
we divide the sky into regions where we assume that the dust- 
H i correlation is constant. For the sake of convenience, we use 
regions of approximate size 13 (52) deg 2 defined by the Healpix 
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pixels at resolution N s iA e = 16 (8) that are outside the Galactic 
mask. We test that our conclusions do not depend on this resolu- 
tion. 

The details of our procedure is as follows. We subtract 
the 143 GHz Planck temperature map from each of the 217- 
857 GHz temperature maps to remove the CMB signal (this 
CMB subtraction is only done for the purposes of creating the 
dust template). We upgrade each of the = 512 LAB maps 
compiled in Land & Slosar (2007) to the Planck map resolu- 
tion of A^ s ide = 2048. Within each region we then simultane- 
ously fit for the amplitude of each Hi velocity component in 
the CMB-subtracted maps, and use the two coefficients per re- 
gion to assemble a full-sky (minus the mask) dust template for 
each of the 217-857 GHz channels. We smooth each template 
with a Gaussian of FWHM 10' to remove the discontinuity at 
the patch boundaries, and then subtract the template from the 
original (CMB-unsubtracted) Planck maps. 

We note that the accuracy of this procedure would be diffi- 
cult to evaluate for all possible uses of the map, i.e., whether it 
might constitute a robust component separation method remains 
to be demonstrated. However, in the case of our cross-correlation 
analysis the dust-removal requirements are less severe, since the 
dust emission only contributes to our measurement as noise. We 
will describe later in Sect. 5.2 the effect on the cross-spectrum 
of removing this emission, and will place limits on the residual 
Galactic contamination in Sect. 5.3.5. 

3.3. Measuring cross-correlations 

To estimate the cross-correlation between the lensing potential 
and a tracer t, we calculate 

m 

As the CIB has an approximately f dependence and the lensing 
potential has an l" 1 dependence, we multiply the cross-spectra 
by ( 3 , and then bin it in 15 linearly spaced bins between 100 < 
i < 2000. As we will discuss in Sect. 5, modes with C < 100 
are not considered, due to possible lens mean-field systematic 
effects, and modes with i > 2000 are removed due to possible 
extragalactic foreground contamination. We have tested that our 
results are robust to an increase or decrease in the number of 
{-bins. 

We expect the error bars to be correlated across bins to some 
extent, due to pseudo-Q mixing induced by the mask, and be- 
tween frequencies, because the lens reconstruction noise is com- 
mon. In addition, any foregrounds that are present in multiple 
channels will introduce correlated noise. The foreground mask 
will also induce a coupling between different modes of the un- 
masked map. This extra coupling can be calculated explicitly 
using the mixing matrix formalism introduced in Hivon et al. 
(2002). Using this formalism and our best-fit models we have 
evaluated the correlation between different bins of the cross- 
correlation signal for our nominal binning scheme. We find that 
the mask-induced correlation is less than 2 % across all bins at 
all frequencies. We will thus neglect it in our analysis. For this 
reason, and based on the results we obtain from simulations, 
we do not attempt to "deconvolve" the mask from the cross- 
spectrum (see e.g., Hivon et al. 2002) and instead correct for the 
power lost through masking the maps by a single sky fraction, 
/sky, ignoring the mode coupling. 

As will be discussed later in Sect. 6.1, when we fit models 
to the cross-spectrum we will assume that the noise correlation 
between bins can be neglected and that the band-powers are flat. 



3.4. Simulating cross-correlations 

In order to validate our measurement pipeline and to confirm that 
our estimate of the cross-spectrum is unbiased we create simu- 
lated maps of the lensed CMB and CIB that have the expected 
statistical properties. 

Using the Planck only favored ACDM cosmology as de- 
scribed in Planck Collaboration XVI (2013) we generate a the- 
oretical prediction of the lensing potential spectrum using CAMB 
(Lewis et al. 2000), from which we generate 300 maps of <p that 
are used to lens 300 CMB realizations using the approach de- 
scribed in Planck Collaboration XVII (2013). We then use the 
PER best-fit CIB model to generate CIB auto- and CIB-0 cross- 
spectra, from which we create CIB realizations that are correctly 
correlated with (f> in each HFI band. The PER model that we 
use describes the CIB clustering at HFI frequencies using a halo 
approach, and simultaneously reproduces known number count 
and luminosity function measurements. At each frequency we 
add a lensed CMB realization to each of the CIB realizations 
and then smooth the maps using a symmetric beam with the 
same FWHM as the beam described in Sect. 2.1. Once this set 
of realizations has been generated we apply the reconstruction 
procedure described above to produce an estimate of the lens- 
ing potential map, and then calculate the cross-power spectrum 
using our measurement pipeline. 

These simulations will miss some complexities inherent in 
the Planck mission. They do not take into account inhomoge- 
neous and correlated noise, and we do not simulate asymmetric 
beam effects. In addition, we do not simulate any foreground 
components, and we instead take a different approach to deter- 
mine their contribution. While simplistic, we believe that our 
simulations are good enough for the purposes of this particular 
measurement. In Sect. 5 we will discuss possible limitations, as 
well as how we test for systematic effects that are not included 
in the simulations. 

We use the simulated maps to check that our pipeline cor- 
rectly recovers the cross-spectrum that was used to generate the 
simulations. For the ^-bins used in our analysis, we find that the 
recovered spectrum is unbiased (to within the precision achiev- 
able with 300 simulations), and with a noise level consistent with 
expectations. The noise in the recovered spectrum is discussed 
in Sect. 5.1. 

4. A strong signal using Planck HFI data 

We now describe the result of applying our pipeline to our nom- 
inal data set, i.e., the lens reconstruction at 143 GHz and the 
foreground reduced Planck HFI temperature maps with a 40 % 
Galactic mask, which when combined with the point source 
mask and H i mask leaves 30.4 % of the sky unmasked. The re- 
sults are presented in Fig. 3. The error bars correspond to the 
naive scatter measured within each bin. The thin black line corre- 
sponds to the expected CIB-lensing correlation predicted using 
the PER CIB model (the HOD parameters of the PER 217 GHz 
best-fit model were used at 100 and 143 GHz since no CIB clus- 
tering measurement at these frequencies is available). As can be 
seen from these plots, the noise is strongly correlated across fre- 
quencies, especially at the lowest frequencies where the CMB 
dominates the error budget. A detailed analysis of the uncertain- 
ties and potential systematic errors attached to this result is pre- 
sented in Sect. 5. 

As clearly visible in Fig. 3, a strong signal is detected. To 
set a reference point and naively quantify its statistical signif- 
icance when taken at face value, we define a detection signifi- 
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Fig. 3. Angular cross-spectra between the reconstructed lensing map and the temperature map at the six HFI frequencies. The error 
bars correspond to the scatter within each band. The solid line is the expected result based on the PER model and is not a fit to 
these data (see Fig. 16 for an adjusted model), although it is already a satisfying model. In each panel we also show the correlation 
between the lens reconstruction at 143 GHz and the 143 GHz temperature map in grey. This is a simple illustration of the frequency 
scaling of our measured signal and also the strength of our signal as compared to possible intra-frequency systematic errors. 



cance as follows. We count the number of standard deviations as 
the quadrature sum of the significance in the different multipole 
bins: 



s 

1=1 



c 



T0 \ 



AC 



(2) 



For our nominal parameters this gives us 3.6 cr, 4.3 cr, 8.3 cr, 
31 cr, 42 cr, and 32 cr, at, respectively, 100, 143, 217, 343, 545 
and 857 GHz. Note that these numbers include an additional 
20 % contribution to the statistical error to account for mode cor- 
relations (which we discuss in Sect. 5.1), but do not include sys- 
tematic errors or our point source correction. As a comparison, in 
each panel we plot the correlation between the lens reconstruc- 
tion at 143 GHz and the 143 GHz map in grey. This shows the 
frequency scaling of our measured signal and also the strength 
of the signal, as compared to possible intra-frequency systematic 
effects. This will be studied in depth in Sect. 5. 

This first pass on our raw data demonstrates a strong detec- 
tion that is in good agreement with the expected CIB-lensing 
signal. To get a better intuition for this detection, we show in 
Fig. 4 the real-space correlation between the observed tempera- 



ture and the lens deflection angles. This figure allows us to vi- 
sualize the correlation between the CIB and the CMB lensing 
deflection angles for the first time. These images were generated 
using the following stacking technique. We first mask the 545 
and 857 GHz temperature maps with our combined mask that 
includes the 20 % Galaxy mask, and identify 20,000 local max- 
ima and minima in these maps. We also select 20,000 random 
locations outside the masked region to use in a null test. We then 
band pass filter the lens map between £ = 400-600 to remove 
scales larger than our stacked map as well as small-scale noise. 
We stack a 1 deg 2 region around each point in both the filtered 
temperature map and lensing potential map, to generate stacked 
CIB and stacked lensing potential images. We take the gradient 
of the stacked lensing potential to calculate the deflection angles, 
which we display in Fig. 4 as arrows. The result of the stack- 
ing over the maxima, minima and random points is displayed 
from left to right in Fig. 4. The strong correlation seen already 
in the cross-power spectrum is clearly visible in both the 545 and 
857 GHz extrema, while the stacking on random locations leads 
to a lensing signal consistent with noise. From simulations, we 
expect a small off-set (^ 1") in the deflection field. This offset 
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Fig. 4. Temperature maps of size 1 deg 2 at 545 and 857 GHz stacked on the 20,000 brightest peaks (left column), troughs (centre 
column) and random map locations (right column). The stacked (averaged) temperature maps is in K. The arrows indicate the 
lensing deflection angle deduced from the gradient of the band-pass filtered lensing potential map stacked on the same peaks. The 
longest arrow corresponds to a deflection of 6.3", which is only a fraction of the total deflection angle because of our filtering. This 
stacking allows us to visualize in real space the lensing of the CMB by the galaxies that generate the CIB. A small and expected 
offset (-1') was corrected by hand when displaying the deflection field. 



was corrected for in this plot. We have verified in simulations 
that this is due to noise in the stacked lensing potential map that 
shifts the peak. As expected, we see that the temperature max- 
ima of the CIB, which contain a larger than average number of 
galaxies, deflect light inward, i.e., they correspond to gravita- 
tional potential wells, while temperature minima trace regions 
with fewer galaxies and deflect light outward, i.e., they corre- 
spond to gravitational potential hills. 



5. Statistical and systematic error budget 

The first pass of our pipeline suggests a strong correlation of 
the CIB with the CMB lensing potential. We now turn to in- 
vestigate the strength and the origin of this signal. We will first 
discuss the different contributions to the statistical error budget 
in Sect. 5.1, and then possible systematic effects in Sect. 5.2. 
Although the most straightforward interpretation of the signal is 
that it arises from dusty star-forming galaxies tracing the large- 
scale mass distribution, in Sect. 5.3 we consider other potential 
astrophysical origins for the observed correlation. 



5.1. Statistical error budget 

In this section we discuss any noise contribution that does not 
lead to a bias in our measurement. The prescription adopted 
throughout this paper is to obtain the error estimates from the 
naive Gaussian analytical error bars calculated using the mea- 
sured auto-spectra of the CIB and lensing potential. We find that 
these errors are approximately equal to 1 .2 times the naive scat- 
ter within an ^-bin, and we will sometimes use this prescription 
where appropriate for convenience (as will be stated in the text). 
This is justified in Appendix A where we consider six different 
methods of quantifying the statistical errors using both simula- 
tions and data. The Gaussian analytical errors, AC^*, are calcu- 
lated using the naive prescription 

/ sky (2( +l)A( (ACf f = Cj T Cf + (C j /f , (3) 

where as before /s^y is the fraction of the sky that is unmasked, 
Af = 126 for our 15 linear bins between I = 100 and t = 2000, 
Cj T and Ct* are the spectra measured using the data, and 
is the model cross spectrum. This last term provides a negligi- 
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Fig. 5. Naive analytical estimates of the contribution to the cj** 
variance as a function of multipole and frequency as given in 
Eq. 4. We assume the same bin sizes as in Fig. 3. The different 
lines are the contribution to the analytical error from the signal 

only: C^Cf m +(cf B *f (green), noise only: Cf N Cf B ' N (blue), 

and the mixed signal and noise terms: C^C^ IB ' N (yellow) and 

C^' N C^ IB (orange). The total contribution is the solid black line, 

and the theory spectrum, (C^ IB ^ , is the dashed line. 



ble contribution due to the large noise bias on Cy, as we now 
describe. 

The statistical error has two sources: instrumental and astro- 
physical. The measured auto-spectra in Eq. 3 contain a signal 
and noise contribution: C; 



C xx + C XX ' N . It is informative to 



split the right hand side of Eq. 3 into four pieces: 



rhs — 



c*;cf m + (c c ( j 



-,ciB0y 



+ £#,N£,CIB,N 



W^CIB.N 



+ cyc e 



^.N^CIB 



(4) 



+ cr c 



Here the first term is a signal-only piece, the second is a noise- 
only piece, and the remaining two terms are mixed signal and 
noise pieces. To discuss the relative importance of these terms, 
we will use for the signal terms the model spectra, and for the 
noise terms we subtract the model spectra from the measured 
spectra: C XXN = Cf x - Cf x . With this definition, the noise will 
include the instrumental contribution, as well as other astrophys- 
ical signals including the CMB, which we do not remove from 
our data for reasons previously mentioned. We show the differ- 
ent terms in Fig. 5. Up to 353 GHz the measured temperature 
spectrum, Cj 1 , is dominated by the CMB at low t and the instru- 
mental noise at high i. At higher frequencies Galactic emission 
dominates at low I and the CIB at high I. For all frequencies up 
to 353 GHz the dominant contribution to the errors in our signal 
comes from the noise-only term (in blue), which is proportional 
to the temperature noise spectrum. At 353 GHz and above the 
mixed signal-noise term CT Cf 1B (orange) becomes important 
and is the largest contribution at 545 and 857 GHz at high I. 

5.2. Instrumental and observational systematic effects 

A number of systematic errors affect the Planck HFI analysis 
and we briefly discuss some of them here. A more complete dis- 
cussion can be found in Planck HFI Core Team (2011b). We 



will illustrate how the very nature of our measurement, a 3 -point 
function, makes it particularly robust to many systematic effects, 
and we will check for their signatures using null tests. For exam- 
ple, there is no noise bias in the 3-point measurement, and many 
effects that can lead to biases in the auto-spectrum of <p do not 
affect us. 



5.2.1 . Potential sources of systematic error 

We begin by describing our knowledge of known systematic ef- 
fects, before discussing others that could bias our result. To ac- 
count for an error in the calibration of the temperature maps, we 
simply add in quadrature a calibration uncertainty to our error 
bars. In Sect. 5.2.2 we use null tests to check that these errors 
are consistent with the data. In addition we use the null tests to 
search for evidence that the calibration has changed between sur- 
veys, for example due to gain drifts. We account for beam errors 
in a conservative manner by using a constant error at each fre- 
quency equal to the maximum error in the beam multipoles, B{, 
at any £ (see PIR for details). The Be uncertainty is larger at high- 
l with, for example, values at € = 1500 of 79.5 %, 8.2 %, 0.53 %, 
0.95 %, 0.31 %, and 0.70%, at 100-857 GHz, respectively. The 
calibration error is therefore larger than the beam error at all 
( between 217 and 857 GHz but smaller at high t in the 100 
and 143 GHz channels. We add the beam error in quadrature in 
an I- and frequency-dependent manner. As discussed in Planck 
Collaboration XVII (2013), uncertainties in the beam transfer 
(as well as the fiducial CMB power spectrum CjT) propagate 
directly to a normalization uncertainty in the lens reconstruc- 
tion. Based on the beam eigenmodes of Planck Collaboration 
VII (2013), it is estimated in Planck Collaboration XVII (2013) 
that beam uncertainty leads to an effective normalization uncer- 
tainty of approximately 0.2% and 143 and 217 GHz, and 0.8% at 
100 GHz. To be conservative, on top of the calibration and beam 
error we will add in quadrature a 2 % uncertainty on the overall 
lens normalization. 

CMB lens reconstruction recovers modes of the lensing 
potential through their anisotropic distorting effect on small- 
scale hot and cold spots in the CMB. The quadratic estima- 
tor, which we use to reconstruct the lensing potential is opti- 
mized to measure the lensing induced statistical anisotropy in 
CMB maps. However, other sources of statistical anisotropy, 
such as the sky mask, inhomogeneous noise, and beam asymme- 
tries, produce signals, which can potentially overlap with lens- 
ing. These introduce a "mean-field" bias, which we estimate us- 
ing Monte Carlo simulations and subtract from our lensing esti- 
mates. Innaccuracies in the simulation procedure will lead to er- 
rors in this correction, particularly if the correction is large. The 
mean-field introduced by the application of a Galaxy and point- 
source mask, for example, which can be several orders above 
magnitude larger than the lensing signal at £ < 100. This is 
discussed further in Appendix B of Planck Collaboration XVII 
(2013). The mask mean-field is a particular concern for us be- 
cause it has the same phases as the harmonic transform of the 
mask. If our masked CIB maps have a non-zero monopole, for 
example, it will correlate strongly with any error in the mask 
mean-field correction. For this reason we do not use any data 
below { — 100 in our analysis. 

To summarize, we do not expect these known systematic ef- 
fects to be present at a significant level. Nevertheless, we still 
perform a set of consistency tests that would be sensitive to them 
or other unknown effects. 
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5.2.2. Null tests 

The Planck scanning strategy, its multiple frequency bands and 
its numerous detectors per frequency, offer many opportunities 
to test the consistency of our signal (see Sect. 2.1). We focus on 
such tests in this section. Our aim is to reveal any systematic ef- 
fects that could lead to a spurious correlation. For all of the tests 
presented, we will quote a^ 2 value as well as the number of de- 
grees of freedom (A^of) as a measure of the consistency with the 
expected null result. Throughout this section, black error bars 
in plots will correspond to the measured scatter within an i bin 
multiplied by 1.2, as was justified in Sect. 5.1 and Appendix A, 
and will also include a CIB calibration error and a beam error, 
while the coloured boxes correspond to the analytical errors of 
the corresponding signal (i.e., not the difference corresponding 
to the null test). Plotting these two error bars illustrates how im- 
portant any deviation could be to our signal. Throughout this 
section, we will illustrate our findings with the 545 GHz corre- 
lation, since it is our prime band for this measurement, but our 
conclusions hold at other frequencies. 

The first test we conduct is to take the temperature differ- 
ence between the two half-ring (HR) maps to cancel any sig- 
nal contribution, and therefore investigate the consistency of our 
measurements with our statistical errors on all time scales. We 
null the temperature maps and correlate with our nominal lens- 
ing map. The results are shown in the left panel of Fig. 6. We see 
a significant deviation from null only when considering survey 
differences. This particular failure can probably be explained by 
apparent gain drifts due to nonlinearity in the analog-digital con- 
version (Planck Collaboration VI 2013; Planck Collaboration 
VIII 2013), not yet corrected at this frequency. Note however 
that the predicted variation is about 1 % while the deviation from 
null would call for a variation of 1.5-2%. But in any case, its 
amplitude is too small to significantly affect our measurement. 

We see a significant deviation from null only when con- 
sidering survey differences. This particular failure can proba- 
bly be explained by apparent gain drifts due to nonlinearity in 
the analog-digital conversion (Planck Collaboration VI 2013; 
Planck Collaboration VIII 2013), not yet corrected at these fre- 
quency. But in any case, its amplitude is too small to significantly 
affect our measurment. 

The second test uses multiple detectors at a given frequency 
that occupy different parts of the focal plane. These detector sets 
are used to construct the "detset" maps that were described in 
Sect. 2. 1 . The two "detset" maps are subtracted and then corre- 
lated with our nominal lens reconstruction. This test is particu- 
larly sensitive to long term noise properties or gain variations, as 
we do not expect these to be correlated from detector to detector. 
Since this detector division breaks the focal plane symmetry, it 
is also a good check for beam asymmetry effects. Here again, we 
do not find any significant deviation, as illustrated in the middle 
panel of Fig. 6. 

The third test we conduct makes use of the redundant sky 
coverage, using multiple surveys to cancel the signal. As above, 
we null the temperature signal and correlate with the nominal 
lens reconstruction. This test is particularly sensitive to any long 
term, i.e., month timescale drifts that could affect our measure- 
ment. It is also a good test for any beam asymmetry effects, as 
individual pixels are observed with a different set of orientations 
in each survey. Since only the first two surveys are complete 
for this particular data release, we only use the two full survey 
maps to avoid complications with the partially completed third 
survey. Here again, we do not find any significant deviation, as 
illustrated in the right panel of Fig. 6. 
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Fig. 6. Null tests at 545 GHz. Left: difference spectra obtained by 
nulling the signal in the HR temperature map before correlating 
it with our nominal <p reconstruction. Middle: temperature signal 
nulled using different detectors at 545 GHz. Right: temperature 
signal nulled using the first and second survey maps. The black 
error bars correspond to the scatter measured within an ^-bin, 
while the coloured bands correspond to the analytical estimate. 
Except for the survey null test (see text for details), these tests 
are passed satisfactorily except, as illustrated by the quoted x 2 
and iVdof . thus strengthening confidence in our signal. 
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Fig. 7. Left: difference between the cross-spectra measured us- 
ing the 20 % Galactic mask (20 % is the unmasked sky frac- 
tion) from that measured with our default 40 % Galactic mask. 
Middle: spectra obtained when differencing the 60 % and 40 % 
Galaxy mask measurements. For both left and middle panels and 
all Galactic masks, the same point source and H i masks are used, 
which removes an additional fraction of the sky. Right: differ- 
ence between the cross-spectra calculated with the H i cleaned 
temperature maps from those with no Hi cleaning. This cross- 
spectrum is thus the correlation between the H i template and the 
<p reconstruction. The error bars are calculated in the same way 
as in Fig. 6. Again, the null tests are passed with an acceptable 

x 2 - 



To conclude, this first set of stringent consistency tests have 
shown that there is no obvious contamination of our measure- 
ments due to instrumental effects. In addition, the reasonable 
X 2 '/Ndaf obtained gives us confidence in our statistical noise eval- 
uation. Although we measure the noise directly from the data, 
this success was not guaranteed. 



5.3. Astrophysical contamination 

We now turn to possible astrophysical biases to our measure- 
ment. We will discuss successively known astrophysical contam- 
inants that can either come from Galactic or extragalactic origin. 
Once again, besides our knowledge of these signals, we will rely 
heavily on consistency tests made possible by having multiple 
full sky frequency maps. 
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Fig. 8. Left: difference between cross-spectra calculated using 
the lens reconstruction at 100 GHz with the nominal 143 GHz 
reconstruction. We see an overall shift, which leads to a high 
reduced^ 2 . This shift can be explained by the expected overall 
normalization uncertainties of the 100 GHz and 143 GHz recon- 
structions. While this uncertainty is not included in the x 2 or 
the solid bars, it is included later in our analysis in Sect. 6.1. 
Middle: same as the left panel, but the 217 GHz reconstruction 
is used instead of the 100 GHz reconstruction. Right: difference 
between cross-spectra when we consider the 143 GHz lens re- 
construction calculated with a less restrictive Galaxy mask (that 
excludes only 20 % of the sky) and the nominal reconstruction 
mask that excludes 40 % of the sky. 



5.3.1. Galactic foregrounds 

Galactic foregrounds have two possible effects on our measure- 
ment. The first is the introduction of an extra source of noise. The 
second is that contamination of the lensing reconstruction by any 
Galactic signal, e.g., synchrotron, free-free or dust, which could 
then correlate with foreground emission present in the tempera- 
ture maps, remains a source of bias that has to be investigated. 
We will show that this bias is small. To do so, we take three ap- 
proaches. We first investigate various Galactic masks, then per- 
form the lensing reconstruction at various frequencies, and fi- 
nally investigate the effect of a dust-cleaning procedure. 

First, we consider two additional masks, either more aggres- 
sive or more conservative than our fiducial one. Both were intro- 
duced in Sect. 2.1. The first one leaves approximately 20% of 
the sky unmasked, while the second one leaves approximately 
60 % of the sky. Given the strong dependence of Galactic fore- 
grounds on Galactic latitude, any Galactic contamination should 
vary strongly when we switch between masks. Comparing the 
measurements using these masks with our fiducial 40 % mask in 
the left and centre panels of Fig. 7, we do not see any substantial 
deviation from our fiducial measurements. This excludes strong 
Galactic contamination of our results. 

Second, we perform the lens reconstruction at 100 and 
217 GHz, different from the fiducial frequency of 143 GHz, and 
compare their correlation with the temperature maps. Given the 
strong dependence of the Galactic emission with frequency, T oc 



detailed in Sect. 3.2 that aims to remove the H i-correlated dust 
component. This procedure leads to a decrease in the variance 
measured outside the mask of 22, 65, 73 and 73 % in the 217, 
353, 545 and 857 GHz maps, respectively. This frequency de- 
pendence is expected given the dust scaling. However, in Fig. 7, 
where we show the differences between the cleaned and non- 
cleaned cross-spectra, we observe that the large scale H i clean- 
ing, even though it makes a substantial impact on the power 
within our map, only makes a small change at low-^ in the cross- 
spectrum, as well as reducing the noise at all multipoles. If we 
quantify the effect of our "local" Hi cleaning on the detection 
significance level computed using only statistical errors, we find 
that the significance is increased by 4, 4, 28, and 36 % at 217, 
353, 545 and 857 GHz, respectively. Also, not surprisingly, we 
observe that for frequencies up to 353 GHz where the statistical 
errors are dominated by the CMB, the H i cleaning has almost no 
effect on the cross-spectra. From the three studies in this section 
we conclude that there are no obvious signs of Galactic fore- 
ground contamination in our cross-correlation. 

5.3.2. Point source contamination 

We now discuss another well-known potential source of contam- 
ination, namely the contribution of unresolved point sources vis- 
ible either through their radio or dust emission. Our concern is 
that a correlation between a spurious lens reconstruction caused 
by unresolved point sources can correlate with sources in the 
temperature map. Although in Sect. 5.3.1 our null test using lens 
reconstructions at different frequencies suggests that unresolved 
point sources are not an obvious contaminant, we will now per- 
form a more detailed test designed specifically to search for point 
source contamination. Following Smith et al. (2007); Osborne 
et al. (2013), we will construct a point source estimator designed 
to be more sensitive than the lensing estimator to point source 
contamination. Our focus here will be on possible contamina- 
tion from the point source shot-noise bispectrum. In Sect. 5.3.5 
we will discuss contamination from a scale dependent bispec- 
trum. 

Our (unnormalized) quadratic estimator, which is designed 
to detect point source contributions is given by 
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where is the inverse-variance filtered sky map. This estima- 
tor is simply the square of the inverse-variance filtered sky map, 
which is a more sensitive probe of point sources than the stan- 
dard lensing estimator. 



In Fig. 9 we plot the cross-spectrum of (0 (n)) LM mea- 
sured at 143 GHz and ® V LM for the full set of HFI channels. This 



for synchrotron and T oc v 2 for dust in this frequency range, cross -spectrum is probing the same point source contributions 



any contamination of our signal would have a strong frequency 
dependence. The comparison with the 100GHz (217 GHz) re- 
construction is presented in the left (centre) panel of Fig. 8. The 
right panel shows the difference of the cross-spectra calculated 
using the 143 GHz reconstruction with a more aggressive Galaxy 
mask (20 % instead of 40%), to reduce possible Galactic contam- 
inants in the reconstruction, and the nominal reconstruction. The 
three differences are consistent with null as demonstrated by the 
quoted^ 2 and iV dof . 

Third, we investigate more specifically how cirrus, the domi- 
nant Galactic contaminant for our higher frequency channels, af- 
fects our measurements. We rely on the dust cleaning procedure 



that could bias our estimates of cJ ', however with a greater 
signal-to-noise ratio. 

There is one complication here, which is that just as lens 
reconstruction may be biased by point source contributions, the 
point source estimator is correspondingly biased by lensing. The 
bias to the plotted cross-spectra is given by 
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J dhY[ im (n)Ye 2m2 (h)Ye 3m (h), C e is the un- 



where - j u,^ (im 

lensed CMB spectrum, Cf is the spectrum of @ 143 (n), /^"f 3 
is as defined in Okamoto & Hu (2003), and we have used the 
fact that for our inverse-variance filtering, &c„ t w &c m IC l ° l . We 

have calculated this contribution using our measured Cjj and 
subtracted it from the data points of Fig. 9. 

We can consider the effect of shot noise on this cross- 
spectrum. With the shot-noise bispectrum defined by 



(s>) 



(71 



the bias to the plotted cross-spectrum is given as 



((® H \h)) 2 LM ®l M ) s , 



(8) 

This bias is plotted for best-fit values of (S 3 ) as the black lines in 
Fig. 9. To minimize systematic effects that might bias the (S 3 J 
estimator, we have estimated S 3 from the spectra of Fig. 9 be- 
tween multipoles between i = 500 and 2000. The fitted ^S 3 ) 
amplitudes are given in Table 1 . 

These amplitudes match our expectations, for example 
see Planck Collaboration XIII (201 1). We observe a decrease in 
the amplitude of the point source contribution going from 100 to 
217 GHz, which corresponds to a dominant contribution from ra- 
dio point sources. We do not see any evidence of a dusty galaxy 
contribution to the shot-noise bias. These conclusions have been 
verified using less restrictive point source masks that cover fewer 
sources. 

With estimates of S 3 in hand, we estimate a corresponding 
bias to C^, given by 



/ j / j £tot£tot Lr (ie 2 L W 

i\m\ i^nii t\ £2 



where W^[ M = (- 1 y*(£™ Q, + 1'^ 

ft^ defined in Planck Collaboration XVII (2013). We show this 
contribution later in Fig. 1 1 as the dotted line. While non-zero, 
we see that the point source shot noise contribution is always 
negligible in the t range we consider, except at lower frequencies 
where the radio point sources are important (but still not strong 
enough to lead to any clear signal in the cross -spectra). 
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5.3.3. SZ contamination 

A fraction of CMB photons travelling from the surface of recom- 
bination are scattered by hot electrons in galaxy clusters. In the 



Table 1. Point source estimator. The measured quantity (S 3 ^, as 
defined in Eq. 9 is given as a function of frequency. 



Frequency 


(53) 


(No. ofcr) 


[GHz] 


[xW 9 fiK 3 ] 




100 


11.7 ± 5.8 


(2.0) 


143 


4.3+ 1.8 


(2.3) 


217 


3.7 ± 1.6 


(2.2) 


353 


6.1 ± 3.9 


(1.6) 


545 


-79 ± 39 


(-2.0) 


857 


(-1.9 ± 2.1) x 10 3 


(-0.9) 
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Fig. 10. Frequency spectrum of our cross-spectra averaged 
within {-bins (black points with associated error bars). The light 
shaded regions correspond to the HFI frequency bands. The solid 
black curve corresponds to the best-fit CIB assuming a Gispert 
et al. (2000) spectrum, while the dot-dashed line assumes a 
Fixsen et al. (1998) spectrum. The dashed black line corresponds 
to the best-fit model when allowing for an SZ component in ad- 
dition to the Gispert et al. (2000) CIB shape. The blue dots cor- 
respond to the associated absolute value of the best-fit SZ com- 
ponent. We conclude from this plot that the SZ effect is not an 
important contaminant. 



most massive clusters approximately 1 % of CMB photons pass- 
ing through them get scattered. On average, their energy will be 
increased, which leads to a measurable spectral distortion. This 
is the so-called thermal Sunyaev-Zeldovich (SZ) effect (Sunyaev 
& Zeldovich 1970). At the location of a galaxy cluster the CMB 
appears colder at frequencies below about 220 GHz and hotter 
at higher frequencies, with a temperature change proportional 
to the cluster optical depth to Compton scattering and to the 
electron temperature. Since hot electrons in clusters also trace 
the large scale matter potential that is traced by CMB lensing, 
we expect an SZ-induced contamination in our measurement at 
some level. We will show below that the level of contamination 
is negligible. In these calculations we ignore the small relativis- 
tic corrections to the thermal SZ spectrum (e.g., Nozawa et al. 
2000). We also ignore the kinetic SZ signal coming from the bulk 
motion of hot electrons in clusters, since it is subdominant to 
the thermal signal (Sunyaev & Zeldovich 1980; Reichardt et al. 
2012; Hand et al. 2012). 

The frequency dependence of the SZ signal in our map de- 
pends on the detector bandpasses and is 



/(v) = 



J dv h(v) g(v) 
J dv h(v) ' 



(10) 



where h(v) is the detector bandpass and g(v) is the SZ fre- 
quency dependence, which in the non-relativistic limit is g(v) = 
x (e x + l)/(e* - 1) - 4, with x = hv/k B T C MB- The effect of the 
bandpass only makes a large difference at 217 GHz near the null 
of the SZ signal. The thermal SZ affects our measurement in two 
ways. First, since the SZ emission in our maps is not a Gaussian 
random field (e.g., Wilson et al. 2012) it introduces a spurious 
signal into our lens reconstruction that will correlate with the 
SZ signal in our CIB map. As shown in Osborne et al. (2013), 
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Fig. 9. Results from the point source contamination estimator of Eq. (5). The best-fit cross-spectra associated with shot noise are 
plotted in black. We do not show the best-fit at 545 and 857 GHz since the signal-to-noise ratio is low. The grey line is a prediction 
for the bias from the CMB lensing - infrared correlation, and has been subtracted from the spectra (plotted as black points). We 
see that with the subtraction of the bias from CMB lensing, the measured bispectrum-related spectrum is generally consistent either 
with zero, or with the shape expected for shot noise. 



this is well approximated by a Poisson noise term and is thus al- 
ready addressed by our treatment of point sources in Sect. 5.3.2. 
The spurious lensing signal will also correlate with other com- 
ponents in our map such as the CIB. However, we ignore these 
terms since they will be smaller than those that correlate directly 
with the SZ emission. Additionally, a contribution comes from 
SZ emission in our CIB map that correlates with the lensing po- 
tential itself. The latter is the dominant term and we discuss it in 
this section. 

To measure a contribution from the SZ-lensing correlation 
we attempt to separate the SZ and CIB emission based on their 
differing spectral shapes. We consider all frequencies from 100 
to 857 GHz, but we will illustrate this procedure by considering 
only two t bands: I = 300-450; and 1200-1450. The first is well 
inside the linear regime, while the second receives a more im- 
portant non-linear contribution. However, we have checked that 
if we consider different £-bins we obtain similar conclusions. We 
model the signal within each I band as S((v) = ai t ec(y)+ci2,tf(y), 
where c(v) and f(v) are, respectively, the CIB frequency depen- 
dence (as proposed in Fixsen et al. 1998 or Gispert et al. 2000) 
and the SZ frequency dependence obtained from Eq. 10. For 



each £ band, we will solve for a\j and a%( minimizing the as- 
sociated x 1 while forcing both amplitudes to be positive. As an 
approximation to the error in each multipole band we calculate 
the scatter of the signal within the band and multiply it by 1 .2, 
as discussed in Sect. 5.1. 

In Fig. 10 we show the measured frequency spectrum within 
each { band, along with the best-fit SZ-lensing and CIB-lensing 
spectra. For the CIB-only fit with the Gispert et al. (2000) fre- 
quency dependence we find a relatively poor fit in the lowest 
^-bin, x 2 (dof) = 15.5 (5), but an improved fit in the higher l- 
bin, x 2 (dof) = 4.15 (5). Including the SZ component gives Ay 2 
= 0.52 and 1.34 in the low and high t bins for one extra degree 
of freedom. When we use the Fixsen et al. (1998) frequency de- 
pendence we find an improved fit, with x 2 (dof) = 2.25 (5) and 
5.49 (5) in the low and high-f-bins, respectively. Overall, the im- 
provement in the ^ 2 /dof when including the SZ component does 
not justify inclusion of the SZ component in the model, with 
the poor fit driven by the lowest frequency bands where the CIB 
scaling is rather unconstrained. In fact, our measurements might 
constitute the first constraints to date on this scaling. From these 
results we conclude that including the SZ-lensing correlation in 
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our data does not improve the fit in the I range of interest to us 
and thus we do not consider it necessary to correct for. 

As an extra validation of this result, we now verify its consis- 
tency with current models of the CIB and SZ emission. For this 
purpose, we use the calculation of the correlation from Osborne 
et al. (2013), based on Babich & Pierpaoli (2008), which models 
the SZ emission as a statistically isotropic signal modulated by 
a biased density contrast, where the bias depends on the clus- 
ter mass and redshift. To obtain an estimate of the contribution 
to the cross-spectrum at 217-857 GHz we assume that the mea- 
sured cross-spectrum at 143 GHz is entirely due to thermal SZ 
emission (note that we do this to find what we believe to be an 
upper limit on the SZ contribution at 217-857 GHz; for the rea- 
sons stated above we do not expect the 143 GHz correlation to 
be due to SZ). Since the SZ signal at 143 GHz gives a decrement 
in the CMB, and the CIB emission gives an enhanced signal, it 
is possible that this approach could still underestimate the SZ 
signal. We find that in order to fit the cross-spectrum at 143 GHz 
using only the SZ-lensing correlation requires an amplitude of 
(2.4 + 1.6) times our calculated SZ-lensing cross-spectrum. In 
Fig. 1 1 the dashed line shows the magnitude of this SZ signal 
scaled to each frequency using Eq. 10. The small contribution it 
makes at 217-857 GHz further suggests that we can neglect this 
component. At 217 GHz the signal is negative, while at higher 
frequencies it is positive. 

5.3.4. ISW contamination 

The Integrated Sachs-Wolfe (ISW) effect describes the redshift- 
ing (blueshifting) of photons travelling through gravitational po- 
tential wells (hills) that decay as the photons travel through 
them (Sachs & Wolfe 1967). The induced modulation of the 
CMB mean by the gravitational potential generates CMB fluc- 
tuations that correlate with the lensing potential, which also 
traces out the gravitational potential perturbations (Seljak & 
Zaldarriaga 1999; Goldberg & Spergel 1999; Lewis et al. 201 1). 
Note that because the mean of the CIB is relatively much smaller 
than its fluctuation, the ISW induced CIB fluctuations make a 
negligible change to total CIB anisotropy. The CMB ISW in- 
duced signal has the same frequency dependence as the CMB 
and so is only a significant contaminant for us at low frequen- 
cies. We evaluate this signal using a theoretical calculation per- 
formed in CAMB (Lewis et al. 2000). The results are shown as the 
solid line in Fig. 1 1. It is a negligible contribution at all frequen- 
cies, except in the lowest ^-bin of the lowest frequencies, where 
the measured cross-spectrum is consistent with zero. 

5.3.5. CIB bispectrum 

Having calculated the bias from the point source shot noise in 
Sect. 5.3.2, we now discuss a more complicated form of the 
unresolved point source 3-point function that could be present 
in our data, namely the clustering contribution. While unknown 
(although the first detection was recently reported in Crawford 
et al. (2013)), the CIB bispectrum is potentially a direct contam- 
inant to our measurement. Because of the non-linear clustering 
of DSFGs (PER), it has to exist. But because of the very large 
redshift kernel that characterizes the CIB, this non-Gaussian ef- 
fect will be washed out, reducing its importance. Nevertheless, 
we ought to study it carefully. 

If important, this effect would show up as a departure of the 
data from the best-fit curve in Fig. 9, since the best-fit model 
that we used assumes only a Poissonian shot-noise contribution. 



2000 




500 1000 1500 2000 

I 



Fig. 12. Cross-spectrum of the 545 GHz lens reconstruction 
correlated with the 545 GHz temperature map with different 
Galactic masks. The legend gives the visible sky fractions. The 
solid line represents the analytic unclustered shot-noise contri- 
bution fit to the / sky = 0.09 points above I = 1300. 



We do not see any significant deviation in Fig. 9. Still, in or- 
der to isolate this effect we create cross-spectra with increased 
sensitivity to the clustered point source signal. We do this by 
calculating the lens reconstruction at 100 GHz and 545 GHz, 
where, respectively, the radio and dusty point source contribu- 
tion is stronger. The 545 GHz map has a much larger Galactic 
dust signal than our nominal 143 GHz map. However, unlike in 
our fiducial estimates, here we do not attempt to project out dust 
contamination from the map used to perform our lens recon- 
struction as this would also remove some of the CIB signal in 
the bispectrum. As was found in Sect. 5.3.1 the cross-correlation 
between the 100 GHz reconstruction and the 100 GHz temper- 
ature map does not show any large difference with the cross- 
spectrum obtained using the 143 GHz signal. We are thus not 
sensitive enough to detect a bias from the clustering of radio 
sources using this method. However, we do detect a strong cross- 
correlation between the raw 545 GHz lens reconstruction and 
the 545 GHz temperature map. This cross-spectrum is shown 
in Fig. 12 for three different Galaxy masks. The line shows the 
point source shot-noise template derived in Sect. 5.3.2, fit to the 
cross -spectrum with the 10% Galaxy mask at t above 1300. If 
the signal were entirely due to extragalactic point sources, then 
the signal would be independent of masking, and we do see a 
convergence of the signal at high I as the size of the Galactic 
mask is increased. At low i, however, there is a large Galactic 
contribution and the convergence with the reduced mask size 
is less clear. We thus conclude that a strong contribution from 
Galactic dust is present in this measurement at all £. 

We do not attempt to calculate accurately the shape of the 
clustering contribution to the CIB bispectrum here, since it is 
beyond the scope of this work, even though a simple prescrip- 
tion for it has recently been proposed in Lacasa et al. (2012). 
To separate the Galactic from non-Galactic contributions in our 
bispectrum measurement is difficult, even if a strong Galactic 
signal is clearly present, given the strong dependence of the sig- 
nal on variations of the Galactic mask in Fig. 12. However, the 
combination of dust cleaning that we perform in our nominal 
pipeline, coupled with the fact that our nominal pipeline uses 
the 143 GHz map for lens reconstruction, means that we do not 
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Fig. 11. Foreground components at each frequency. The data points and error bars show our results. The dashed line is an estimated 
upper limit on the magnitude of the SZ contamination derived in Sect. 5.3.3. We show the absolute value of this contribution, which 
is negative at frequencies less than 217 GHz. The dot-dashed line is the extragalactic point source contribution, with an amplitude 
measured from our data as derived in Sect. 5.3.2. Again we show the absolute value, with the signal being negative below t ~ 1200. 
The oscillating solid line corresponds to the calculated ISW contamination. 



observe any dependence with masking in our measurement, as 
seen in Fig. 7. Because of this, the CIB bispectrum is unlikely to 
be a large contribution to our measurement. Furthermore, even 
if we were to assume that all of the signal seen in Fig. 12 was 
extragalactic in nature, using the Gispert et al. (2000) frequency 
scaling for the CIB (also appropriate for the Galactic dust in fact, 
Planck Collaboration XXIV 201 1), the roughly - 1700 yuK.sr ob- 
served at I — 400 for the 40 % Galactic mask would only lead 
to a -0.02ji/K.sr signal in Fig. 11, which is an order of magni- 
tude smaller than our measured signal. To conclude, although 
our analysis does not lead to a clean measurement of the CIB 
bispectrum, we can safely assume that it is not a contaminant to 
our measurement. 



nificant deviation gives us confidence in our error budget. As 
described in Sect. 5.2 we add to them an overall calibration un- 
certainty, beam uncertainty, and lens normalization uncertainty, 
consistent with the Planck data processing paper (Planck HFI 
Core Team 2011b). We gather the measured band-powers in 
Table 2, along with our statistical and systematic errors. These 
band-powers have been corrected for the point source compo- 
nent measured in Sect. 5.3.2, whose amplitude is also given in 
Table 2. 

Once all systematic effects are factored in, we claim a de- 
tection significance of 3.6 (3.5), 4.3 (4.2), 8.3 (7.9), 31 (24), 42 
(19), and 32 (16) <x statistical (statistical and systematic) at 100, 
143, 217, 353, 545 and 857 GHz, respectively. 



5.4. Final statistical and systematic error budget 

Throughout the suite of tests for instrumental and observational 
systematic errors presented in Sect. 5.2, as well as the suite 
of tests for possible astrophysical contaminants presented in 
Sect. 5.3, we have established the robustness of our measure- 
ment. The fact that our consistency tests do not lead to any sig- 



6. Interpretation and discussion 

The correlation we have investigated leads to a very strong sig- 
nal at most HFI frequencies. After a thorough examination of 
possible instrumental and astrophysical origins, we interpret it 
as originating from the spatial correlation between the sources 
of the CIB and the matter responsible for the gravitational de- 
flection of CMB photons. In this section, we build on this result 
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Table 2. Cross-spectrum detection band-powers. All values are in units of fiK.sr. The extragalactic foreground contribution, C f ° K 
has been removed from Q . Both statistical and systematic errors are given (see Sect. 5 for details). 
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and interpret the measurement using both angular and frequency 
information. 

Before doing so, we highlight the spectral information con- 
tained in the signal. It was shown in PER that both the fre- 
quency spectrum of the CIB mean and fluctuation rms are well 
approximated by the two modified blackbody spectra proposed 
by Fixsen et al. (1998) and Gispert et al. (2000), with a slight 
preference for the latter. We expect our measurements to follow 
the same spectral energy distribution (SED). Following the pro- 
cedure outlined in Sect. 5.3.3, we plot in Fig. 10 the best-fit CIB 
component with either a Fixsen et al. (1998) - dot-dashed black 
line - or Gispert et al. (2000) - solid black line - SED. We do 
so for two ^-bins. We can see that, indeed, for a given (-bin, 
our measurements qualitatively follow the expected CIB spec- 
trum. Unlike PER, we do not find a preference for the Gispert 
et al. (2000) shape in our low £-bm and only a slight one in 
the high £ bin. It is worth emphasizing that by carrying out a 
cross-correlation measurement we can obtain constraints at the 
lowest frequency, which is usually heavily contaminated by the 
CMB (see PIR for discussion). This is particularly interesting, 
because these measurements are simultaneously the most sensi- 
tive to high-z star formation processes and the most discrepant 
with either of the SEDs, i.e., they are both systematically low by 
about 0.5 cr. The models presented in this section will allow us 
to use both the spectral dependence and the relative amplitudes 
of the ^-bins that was lost in Fig. 10. We now describe the gen- 



eral methodology we will use, before describing our models in 
detail. 

6.1. Model comparison methodology 

For the purpose of model fitting, we will utilize both the CIB- 
lensing cross-spectra measured in this paper and the CIB auto- 
spectra obtained from PIR. We use the CIB-lensing cross-spectra 
for two purposes: to improve constraints on the model parame- 
ters; and to provide a consistency test of models fit to the CIB 
auto- and frequency cross-spectra alone. As will be seen in PIR, 
the cross-spectra of the CIB at different frequencies provide 
powerful constraints on the CIB emissivity. 
We use the log-likelihood, 

ln£(p) = Z [ C < " & M WT [3 - ' 

11} 

( n ) 

where C and C are the data and theory spectra with parameters p, 
the i and j indices denote the type of spectra (e.g., 100 GHz x <p 
or 100 GHz x 100 GHz), and N is the covariance matrix that in- 
cludes both statistical and systematic errors. We make the ap- 
proximation that the covariance matrix is diagonal, i.e., we treat 
the errors for different bins of each auto- and cross-spectrum 
as being uncorrelated. The small (^ 2 %) mask-induced mode- 
coupling between neighbouring bins supports this approxima- 
tion. However, calibration and beam errors (which are corre- 
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lated between the auto- and cross-spectra at a given frequency), 
as well as the lens normalization error (which is also corre- 
lated across spectra) are not accounted for in this approxima- 
tion. In addition, the lens reconstruction has some sensitivity to 
all modes of the temperature maps, and so different <p modes 
are correlated to some degree. We also neglect the fact that the 
contribution to the error from the CIB signal itself (the orange 
line in Fig. 5) is also substantially correlated from frequency to 
frequency. However, our evaluation using simulations suggests 
that these effects are too small to significantly affect our pro- 
cedure. We thus resort to simply adding the beam, calibration 
and normalization uncertainties in quadrature to the statistical 
errors. The posterior probability distributions of model param- 
eters are determined using now standard Markov Chain Monte 
Carlo techniques (e.g., Knox et al. 2001; Lewis & Bridle 2002). 

6.2. Two modelling approaches 

The strength of the correlation signal should come as no sur- 
prise, given our current knowledge of CMB lensing and the CIB. 
The PER model predicts a high correlation between the CIB and 
the lensing potential. As clearly illustrated in Fig. 1, the broad 
overlap of the redshift distributions of the CIB with the lensing 
kernel peaking at z ~ 2-3 leads to a correlation of 60-80 %. This 
is comparable to our measurements at all of the HFI frequencies, 
as illustrated in Fig. 13, 

In models of the cross-correlation, the underlying properties 
we can probe come from a combination of the lensing potential 
and the characteristics of the DSFGs, in particular their emis- 
sivity and clustering properties. Mostly driven by linear physics, 
the former is well understood theoretically, as confirmed by re- 
cent observations (Smith et al. 2008; Hirata et al. 2008; Das et al. 
2011; van Engelen et al. 2012). Assuming the currently favored 
ACDM cosmology, we can consider it to be known to better than 
10 % in the multipole range of interest to us, an uncertainty dom- 
inated by the uncertainty in the normalization of the primordial 
power spectrum. Given that this is much smaller than the theoret- 
ical uncertainties related to DSFGs, we will fix the cosmology to 
the currently favoured Planck alone flat ACDM model in Planck 
Collaboration XVI (2013), and will focus our analysis on the 
modelling of the DSFGs. 

At a given redshift we model the fluctuations in the mean 
CIB emission, j, as being proportional to the fluctuations in the 
number of galaxies, « g (Haiman & Knox 2000), 



<* j — 

Ha 



(12) 



With this hypothesis, the goal of the CIB modelling becomes 
twofold: first, to better understand the statistical properties of 
the dusty galaxy number density, 6n g ; and second, to reconstruct 
the mean emissivity of the CIB as a function of redshift. 

In this paper we will use two different models of the CIB 
emission. The first model (described in Sect. 6.2.1 and inspired 
by Hall et al. 2010) uses a single bias parameter at all frequen- 
cies with the mean CIB emissivity modelled as a two parameter 
Gaussian. This model is not designed to be physically realistic, 
and furthermore we will marginalize over an arbitrary amplitude 
in this case. Nevertheless, we present this simple model to show 
that our measurements are quite consistent with broad expecta- 
tions of the CIB. The second model, described in Sect. 6.2.2, is 
a natural extension of the Halo Occupation Density (HOD) ap- 
proach used in PER (see also Penin et al. 2011, and references 
therein). But unlike the results obtained in PER we now use a 




2000 



Fig. 13. Cross-correlation coefficients calculated from the model 
<p spectrum and best-fit halo model at each frequency. The CIB 
is a spectacular tracer of CMB lensing, and vice-versa. The data 
points represent the measured cross-correlation divided by the 
best-fit auto power spectra models at 545 GHz. 




Fig. 14. Marginalized 2-D distributions of z c and o~ z for the linear 
bias model, fit to all frequencies simultaneously. The orange dots 
indicate the parameter values at the minimum^ 2 . 



single HOD to describe the spectra at all frequencies. This is 
possible by allowing for deviations from the Bethermin et al. 
(201 1) model (hereafter Bl 1) that was used to fix the emissivity. 

6.2.1. Linear bias model 

As a first pass at interpreting our measurement, we will consider 
a redshift independent linear bias model with a simple paramet- 
ric SED. This model was found to provide a reasonable fit to the 
auto-spectra in the linear regime in PER. Throughout this paper 
we use the Limber approximation, and in this section, since we 
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are using a linear model, we write the relevant angular spectra as 



d X W x 

Jo 



(xW y (x) p ss (k = e/ x ,x), 



(13) 



where X and Y are either the CIB at frequency v or the lensing 
potential 0, the integral is over x, the comoving distance along 
the line of sight, x * 1S the comoving distance to the last scattering 
surface, Pss(k,X) is the matter power spectrum at distance^, and 
the W x functions are the redshift weights for each of the signals 
X: 



W v (x) = b 



ajv(x) , 
X 



(14) 



Here b is the DSFG bias that we assume to be redshift indepen- 
dent, a is the scale factor, j v (x) is the mean CIB emissivity at 
frequency v, as defined in PER, Q,„ is the matter density today 
in critical density units and Hq is the Hubble parameter today. 
We use the Hall et al. (2010) model for the CIB kernel, which 
treats the CIB emissivity as a Gaussian in redshift: 



jv(x) ^ax 2 exp [-(z - z c ) 2 /2<x 2 ] / v(1+z) (15) 

where we use a modified blackbody frequency dependence 

f v{1+z) oc v^BATd). (16) 

We fix the dust temperature to r<j = 34 K, the spectral in- 
dex to (3 — 2 (Hall et al. 2010), and assume a constant bias 
b. We include a single normalization parameter for /', which 
we marginalize over. Since the normalization and bias param- 
eters are degenerate in Eq. 13, if we were to only use the mea- 
sured auto- and cross-spectra this approach would be equivalent 
to marginalizing over a frequency independent bias parameter. 
However, we will further constrain our model using the FIRAS 
data, which breaks this degeneracy. We constrain the z c and <r z 
parameters at each frequency, giving us a total of 13 free param- 
eters. 

For 217-857 GHz, we use the FIRAS measurements of the 
CIB mean intensity from Lagache et al. (2000) as an additional 
constraint to our model. The mean intensity is simply 



-f 

Jo 



d X aj v (x) 



(17) 



Using this equation and the measured FIRAS mean and uncer- 
tainty we calculate a X 2 value and add this to the^ 2 in Eq. 11. 
Since there are no FIRAS constraints at 100 and 143 GHz, as 
well as no CIB auto-spectra measurements, and noisier cross- 
spectra measurements at these frequencies, our constraints for 
the 100 and 143 GHz redshift parameters are weaker than for the 
other parameters. 

The linear bias model considers only linear clustering, and 
so when fitting the auto-spectra we restrict ourselves to t < 500, 
where non-linear contributions are negligible. Because we do not 
consider the high-^ modes, we also neglect the shot-noise con- 
tribution to the auto-spectra. The best-fit model is shown as the 
coloured dashed lines in Fig. 15, with^ 2 values of 13.4, 16.8, 
25.2, 21.8, 9.1, and 9.4 if we break up the X 2 contribution per 
frequency from 100 to 857 GHz, leading to an overall ;f 2 of 95.7 
for A^dof = 59. We see that the model captures some features of 
the data, but we also have evidence it is significantly missing 
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Fig. 15. Measured cross-spectra with the best-fit j reconstruction 
model fit to both the CIB auto- and CIB-lensing cross-spectra 
(solid coloured), and the best-fit linear bias model (dashed 
coloured). The x 2 values quoted in each panel are the contribu- 
tion to the global^- 2 from the data in the panel for the halo model, 
and loosely indicate the goodness of fit (see text for details). 
The one and two-halo contributions are shown as the dashed and 
solid black lines, respectively. A light dashed black horizontal 
line is indicating the zero level. 
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Fig. 16. PIR auto-spectra with the best-fit mean emissivity j re- 
construction model fit for the CIB auto and CIB-lensing cross- 
spectra (solid coloured). The^- 2 values are defined as in Fig. 15. 
The one and two-halo contributions are shown as the dashed and 
solid black lines, respectively, while shot noise is the dot-dashed 
black line. 



some as well. This is perhaps not surprising given the simplic- 
ity of the model. The two-dimensional marginal distributions of 
the Zc, 0"z parameters are shown in Fig. 14. Although we allowed 
for these parameters to be frequency dependent we note that the 
point z c = 1 and cr z = 2.2 is in a region of high probability at all 
frequencies, and gives a redshift distribution for the emissivity 
density roughly consistent with our expectations, rising toward 
Z = 1 due to the x 2 term and then only slowly falling off toward 
higher redshifts. 
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Althought it is useful to see the extent to which such a simple 
model can explain our data, we now turn to make a stronger con- 
nection between the properties of the infrared light and the dis- 
tribution of the underlying dark matter applicable into the non- 
linear regime. 

6.2.2. An extended halo model based analysis 

In this section we use a description of the CIB motivated by the 
halo model, which has been used successfully to describe the 
transition between the linear and non-linear clustering regimes 
for optical galaxies. We use the halo model to attempt to recon- 
struct the CIB emissivity as a function of redshift. This is an ex- 
tension of the approach taken in PER, where the modelled CIB 
emissivity at high redshift was treated as a single bin with the 
amplitude constrained by the data. The goal of this approach is 
to isolate the high-redshift contribution to the CIB, which is dif- 
ficult to probe using observations of individual galaxies, due to 
their low brightness. The power of such an approach is further 
demonstrated in PIR. 

We replace the linear bias used in Sect. 6.2.1 with a halo 
model and an HOD prescription that assigns galaxies to host 
dark matter halos (see PER for references and definitions). It 
allows a consistent description of the linear and non-linear part 
of the galaxy power spectrum and its redshift evolution. Because 
it is built on the clustering of dark matter halos, the halo model 
allows us to describe the clustering of DSFGs and the gravita- 
tional lensing caused by the halos in a consistent way. However, 
it is important to realise that the HOD prescription was devel- 
oped to describe stellar mass within dark matter halos - an ap- 
plication for which it has been thoroughly tested - while here 
we are applying it to star formation within halos. The accuracy 
of this approach needs to be further quantified. However, it pro- 
vides a good phenomenological description of our data as well as 
other CIB measurements, but also of other astrophysical probes 
of the relation between dark matter and light (e.g., Leauthaud 
et al. (2012); Hikage et al. (2012)). 

Unlike the model presented in PER we use a single HOD 
to describe our data at all frequencies. This is made possible by 
allowing for a deviation from the Bll emissivity model. Note 
however that we will still consider the CIB emissivity to depend 
only on redshift and not on the galaxy host halo mass, a simpli- 
fication highlighted in Shang et al. (2012) that will be relaxed 
in the PIR model. The emissivity of the CIB is inhomogeneous, 
due to spatial variations in the number density of galaxies: 

6 j v 6n„ 

Jl(A,z) = — Hn,z) = 6Jn,z). (18) 

jv "g 

Here j(n ,z) is the CIB emissivity at redshift z with mean j(z), 
n g (n,z) is the number density of DSFGs with mean n g (z), 
and 6 g (n,z) is the DSFG overdensity, with power spectrum 

(<5 g (k,z)<5 g (k',z)*) = (27i) 3 6(k-k')P sg (k,z). We calculate this 
power spectrum, including the constituent 1 and 2-halo terms, 
using the procedure described in appendix C of PER, with the 
constraint ff sat = 1, a theoretically favoured value (Tinker & 
Wetzel 2010). We remove the relationship between M sat , a char- 
acteristic satellite mass scale, and M m [ n , the halo mass at which 
a halo has a 50 % probability of containing a central galaxy that 
was enforced in PER (i.e., M sat = 3.3M m ; n ), and allow both M sat 
and M m ; n to vary independently. 

At redshift z < 1 we fix the emissivity to the Bll value, but at 
higher redshift we assume that the emissivity is constant within 
z-bins and solve for the amplitude of the bins. Two factors affect 
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Fig. 17. Marginalized 2-D distribution of log I0 (M m ; n /Mo) and 
log 10 (M sat /M G ) for our overall HOD model when the CIB- 
lensing cross-spectra are combined with the CIB auto-spectra 
and FIRAS measurements. The orange dot shows the best-fit 
value. The contours correspond to 68%, 95% and 99.7% con- 
fidence intervals. 

the number of bins that we choose. The auto-spectra have a j 2 
dependence, and so if the true value of j has a strong z depen- 
dence within a bin then the best-fit emissivity in the bin will be 
difficult to interpret. The best-fit bin values could be significantly 
different from those that would be obtained by binning the true 
emissivity. However, as more bins are used and the number of 
parameters increases, it becomes more difficult to determine the 
best-fit parameters and the parameters will be highly correlated. 
After investigation using simulations, we found that three bins 
was a good compromise, given the expected slow redshift evo- 
lution. The bins are defined by: 1 < z < 1.5; 1.5 < z < 3; and 
3 < z < 7. As in Sect. 6.2.1 we use the FIRAS results at 217- 
857 GHz to add an integral constraint on the emissivity. The CIB 
auto and lensing cross-spectra are (Song et al. 2003): 

cf = fd x w v ( x )w v '(x) p gg (k = ei x ,x)\ 

J r (19) 
c v ; = J dx W v ( X )W^) P Sg (k = llx,x) . 

Since we fix j at z < 1, the model spectra consist of a low red- 
shift part that is independent of the emissivity parameters, and 
a contribution from z > 1 that is proportional to j v j v > for the 
auto-spectra and j v for the lensing cross-spectra. 

Overall, the halo-based model contains two halo param- 
eters that describe the galaxy clustering and are independent 
of frequency, and three j amplitudes at each frequency, giv- 
ing a total of 20 parameters for the six frequencies of interest 
to us. The auto and cross-spectra have a total of 120 ^-bins, 
with four additional FIRAS data points. Solving for the likeli- 
hood described in Sect. 6.1, gives the best-fit models shown in 
Figs. 15 and 16 as solid lines. The^ 2 values in each panel are 
the contribution to the total x 2 from the data within the panel. 
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The combined reduced X 2 is 102.1 for N do{ = 104, indicat- 
ing a good fit. The constraints we find on M sat and M m j n are 
shown in Fig. 17. We force M sat > M mia in the MCMC fit- 
ting procedure, with the dashed line in Fig. 17 showing equal- 
ity. The red cross corresponds to the parameter values that give 
the minimum^ 2 in the fit, and are log 10 (M m i n /M Q ) = 12.18 
and log lo (M sat /M ) = 12.76, which gives M sat /M min = 3.80. 
The mean parameter values are log 10 (M m j n /M ) = 10.53 ± 0.62 
and log 10 (M sat /M G ) = 10.80 + 0.74. The best-fit value of M min 
is consistent with those derived in PER at multiple frequencies, 
even though we now set a sat = 1 and reconstruct the mean emis- 
sivity as a function of redshift. The associated mean emissivity 
parameters are given in Table 3 and displayed in Fig. 18, where 
we also plot the Bll model for reference. As can be seen in 
Fig. 18, we remain consistent with the Bll model in most red- 
shift bins. 



6.3. Interpreting the reconstructed emissivities 

We now illustrate one interesting consequence of this measure- 
ment and show how using the constrained emissivities, j v (z), we 
can estimate the star formation rate (SFR) density at different 
redshifts. Following Penin et al. (201 1), we begin by writing the 
emissivity as an integral over the galaxy flux densities: 



jv(z) = a 



dx 

i— — 

dz 



d 2 N 
' dS v dz 



dS x 



(20) 



Here S v is the flux density, and d 2 N/dS v dz is the number of 
galaxies per flux element and redshift interval. The galaxies con- 
tributing to the CIB can be divided into various populations (la- 
belled as p) based on the galaxy SED, e.g., according to galaxy 
type or dust temperature: 



jv(z) = 



I — 

dz 



d 2 N p 
' dS v dz 



dS % 



(21) 



If we define s v as the flux density of an Lir = L Q source with the 
SED of a given population, i.e., S Y — s v L IR (with L IR in units of 
L Q ), then we can write Eq. 21 as (Penin et al. 201 1): 



jv(z) 
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dz 



dV 
dz 



4R 



d 2 N p 
t/LiR dV 



-dL 



IR- 



(22) 



The contribution to the infrared luminosity density from a given 
population is 



Pir, 



d 2 N p 
dL lR dV 



dLn 



(23) 



We assume a simple conversion between Lir and the star forma- 
tion rate density, psfr, using the Kennicutt constant K (Kennicutt 
1998). Since by definition psfr = KYj P PiR,p, we can rewrite 
Eq. 22 as: 



j v (z) = (1 +z)x 



2 PSFR 



i PlR,p 



UnPlR.p 



(24) 



where the final term in brackets is the effective SED of infrared 
galaxies, which we write as s v , & s- We derive these SEDs follow- 
ing the evolution model of Bethermin et al. (2012) using Magdis 
et al. (2012) templates. The construction of these effective SEDs 
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Fig. 18. Reconstructed mean emissivity, j, for each frequency as 
a function of redshift. The solid line at low z and the dashed line 
at higher z correspond to the Bll model. The Bll emissivity 
model at z > 1 is not used, and is shown only for reference. 
The black error bars correspond to the 68% C.L. while the color 
shading display the full posterior distribution. 
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Fig. 19. Marginalized 1-D distribution of the emissivity in the 
high redshift bin at 353 GHz with (black line) or without (blue 
line) including the CIB-lensing correlation. Its inclusion helps to 
constrain the emissivity at high redshift, transforming an upper 
limit into a detection. 



will be explained in detail in future work. Finally, we obtain the 
conversion factor between mean emissivity and SFR density, 



Psfr(z) = 



K 



(l+z)x 2 (z)S v Mz) 



jv(z) ■ 



(25) 



Using Eq. 25 we find the coefficients for each of the redshift bins 
and frequencies used in Table 3. 

6.4. Discussion and outlook 

In the previous section we described a model that simultaneously 
fits the CIB auto-spectra and the CIB-lensing cross-spectra, at all 
frequencies and with a single HOD prescription. Given that we 
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Table 3. Reconstructed emissivity as a function of redshift and associated star formation rate. At each frequency and for each of 
the three redshift bins the first quantity corresponds to the mean emissivity in the corresponding redshift bin, j(z), in Jy Mpc~' sr _1 , 
while the second corresponds to the SFR density, psfr, in M Q Mpc~ 3 yr~' . 



1 <z< 1.5 



1.5<z<3 



3 <z< 1 



j(z) 



PSFR 



](Z) 



PSFR 



](z) 



PSFR 



100 GHz 
143 GHz 
217 GHz 
353 GHz 
545 GHz 
857 GHz 



7.16±5.77 
12.7±9.60 
11.9+6.33 
116+17.1 
185+106 
193+139 



1.96+1.58 
1.37+0.964 
0.310+0.165 
0.671+0.099 
0.320+0.183 
0.144+0.104 



3.53+3.05 0.655+0.564 5.49+4.78 0.271+0.236 



6.82+5.46 
17.3+7.23 
75.5+27.5 
224+148 
354+212 



0.438= 
0.282= 
0.286= 
0.317= 
0.317= 



0.351 
0.118 
0.104 
0.210 
0.190 



10.5+9.05 
36.6+13.8 
164+47.3 
417+251 
609+359 



0.178+0.153 
0.182+0.068 
0.320+0.092 
0.659+0.396 
1.37+0.809 



50 
40 



30 




I 

Fig. 20. Correlation between the lensing potential and the IRIS 
map at 100 fim using our nominal lens reconstruction. We clearly 
see a correlation and estimate the significance to be 9 cr, ignoring 
possible systematic effects. The solid line represents a simple 
reasonnable prediction for this signal. 



use an emissivity function that is close to the Bll emissivities 
(to within our uncertainties), we expect predictions of the galaxy 
number counts derived from our best-fit emissivity to agree with 
current estimates (Bethermin et al. 2012). The fact that our mea- 
surement is consistent with previous models of the CIB lends 
support to our current understanding of its origin. For exam- 
ple, the characteristic mass scale at which halos host galaxies, 
M m i„, is consistent with values derived in PER, and is consistent 
with, but slightly higher than, the value derived more recently 
in Viero et al. (2012), log 10 (M^/Mg) = 9.9 + 0.5 (although a 
direct comparison could be misleading given the different model 
assumptions). In particular, it is clear that our model has limi- 
tations, some of which have been partially addressed in recent 
work (Shang et al. 2012; Bethermin et al. 2012; De Bernardis & 
Cooray 2012; Bethermin et al. 2012; Viero et al. 2012; Addison 
et al. 2012) and are points of focus in PIR, amongst them the 
mass independence of the emissivity. Another question worth 
further investigation is the dependence of our results on the bin- 
ning scheme chosen for the emissivity, which will be addressed 
in a future paper. 

Given the consistency of our model with the PER results, 
the information added by our cross-spectrum measurement is 
worth quantifying. As an example, we show in Fig. 19 the high- 
est redshift emissivity bin in the 353 GHz band. Adding the CIB- 



lensing cross-spectrum information tightens the constraint on the 
high-redshift part of the emissivity. This statement also holds for 
the other frequencies and stems from the fact that the CMB lens- 
ing kernel peaks at high redshift, making the cross-correlation 
more sensitive to the high-redshift CIB signal than the CIB auto- 
spectrum, as is illustrated in Fig. 1. Although this gain does not 
translate into a substantial improvement in M m ; n , it leads to inter- 
esting constraints on the SFR density, as can be seen in Table 3. 

The results at frequencies above 217 GHz each lead to 
around 2 cr evidence for a non-zero SFR density for 1 .5 < z < 3 
and for 3 < z < 7. The values inferred are consistent with other 
probes of the SFR in these redshift ranges, as compiled for ex- 
ample in Fig. 1 of Hopkins & Beacom (2006). Assuming that 
each frequency is independent, we obtain SFR densities for the 
three redshift bins of 0.423 ± 0.123, 0.292 ±0.138 and 0.226 
± 0.100 M Mpc~ 3 yr _1 , respectively where the errors are 68% 
C.L.. We note that the j distributions are rather non-Gaussian so 
that the 95% C.L. become 0.228, 0.246 and 0.191 respectively. 
This roughly 2cr detection per bin compares very favourably 
with other published measurements. These constraints clearly 
illustrate how this particular correlation can be used to better 
isolate the high redshift component of the CIB and improve our 
constraints on the star formation rate at high redshift. Such con- 
strains will improve with future measurements, in particular if 
we can increase the signal-to-noise ratio in our lower frequency 
channels, where the high redshift contribution is the greatest. 
This will likely require an accurate removal of the CMB, our 
dominant source of noise at low frequencies. A more thorough 
discussion of this possibility will be given in PIR. 

To fully utilize the richness of the correlation will require 
more studies. Future work could involve using more sophisti- 
cated halo models specifically designed to model star formation 
within halos, as well as relaxing some of the assumptions made 
here, such as the mass independent luminosity function. In addi- 
tion the use of map-based methods that enable estimates of the 
galaxy host halo mass by stacking the lensing potential maps is 
worth pursuing, as is the extension to other data-sets. For illus- 
tration purposes, we show in Fig. 20 raw measurements of the 
correlation between our lensing potential map and the IRIS map 
at 100 fim that was introduced in Sect. 2.2.2. We use our nominal 
mask and lens reconstruction, with no H i cleaning performed on 
the IRIS map. We clearly see a strong correlation, whose signif- 
icance we estimate to be 9 cr, ignoring any possible systematic 
effects. To guide the eye we have added a prediction (not a fit) 
based on the HOD model presented in Penin et al. (2011). The 
full analysis of this signal is beyond the scope of this paper, but 
it illustrates possible future uses of the lensing potential map. In 
this case the IRIS wavelength range will help us to isolate the 
low-redshift contribution to the CIB. 
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To conclude, we have presented the first measurement of 
the correlation between lensing of the CMB and the CIB. 
Planck's unprecedented full-sky, multi-frequency, deep survey 
enables us to make an internal measurement of this correlation. 
Measurements with high statistical significance are obtained, 
even after accounting for possible systematic errors. The high 
degree of correlation that is measured (around 80 %) allows for 
unprecedented visualization of lensing of the CMB and holds 
great promise for new CIB and CMB focused science. CMB 
lensing appears promising as a probe of the origin of the CIB, 
while the CIB is now established as an ideal tracer of CMB lens- 
ing. 
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Appendix A: Statistical Errors 

In this section we compare six different methods to estimate our 
statistical errors. This comparison is used to validate our claim 
(presented in the main text) that we can obtain our errors from 
the naive analytical errors calculated from the data, i.e., method 
1 below. The six different methods we compare are: 

1. The naive, Gaussian, analytical errors estimated from the 
data through the measured total power of the T and fields, 
respectively Cj 1 and CZ and the model cross-spectrum C^. 

2. As above but instead of the data maps we use one of our sim- 
ulations of the CIB and CMB temperature maps described in 
Sect. 3.4. The lens reconstruction is obtained from the sim- 
ulated maps using the same procedure that we use for the 
data. 

3. The scatter directly within individual (-bins in the data- 
determined cross-spectrum. 

4. As above but the scatter is measured in each f-bin for each 
simulation realization, and the errors are averaged over 100 
realizations. 

5. The scatter of the bins is calculated using the cross-spectra 
measured from our simulated maps. This is a direct mea- 
surement of the statistical error we require (to the extent that 
our simulations are realistic), and will differ from the scat- 
ter within the ^-bins, for example due to noise correlations 
between different multipoles of the cross-spectrum. 

6. The error in the cross-spectrum of the reconstructed lensing 
potential in simulated maps with the measured temperature 
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Fig.A.l. Ratio of various error estimation procedures to the errors obtained with the data-based analytical estimate. At each fre- 
quency the numerator is given by: (i) the scatter within an ^-bin in simulations (solid black line); (ii) the scatter within an (-bin 
in the data (solid dashed black line); (iii) the scatter of bins across simulated realizations (solid coloured line); (iv) the analytical 
errors calculated from the simulations (dashed coloured line); (v) the scatter across realizations for the cross-correlation between 
the simulated temperature map and the lensing potential reconstructed from the data (coloured dot-dashed line). The grey envelope 
is the precision of the simulated errors expected from 100 simulations (shown as a spread around unity). 



maps. This will only give part of the contribution to the error 
since the temperature maps are fixed, but it is still a useful 
cross-check. 

In Fig. A.l we show a comparison of the errors found from 
our six measurement methods. The precision achievable with 
100 simulations is indicated by the grey envelope. We show 
the errors in each t bin from the different methods divided by 
the data-based analytical estimate. To discuss the implications 
of these results we shall focus on the 100 GHz panel first. The 
scatter measured within an I bin is fairly consistent in the sim- 
ulations (method 4, black solid line) and in the data (method 
3, black dashed line) giving us confidence in our simple simula- 
tion procedure. This rules out important systematic contributions 
and shows that our signal is mostly Gaussian, as expected. Note 
that the consistency with the simulations is not surprising, since 
at low frequencies we are dominated by CMB and instrumen- 
tal noise, which are well understood. In addition, the fact that 
the analytical errors calculated on the simulations (method 2, 
coloured dashed line) are mostly within the grey shaded region, 
and are therefore close to the analytical errors calculated from 
the data (method 1), gives us further confidence in the simula- 
tions. To the extent that the simulations are accurate, the scatter 
of the ^-bins in simulations (solid coloured line) is the error that 
we require. The fact that it is essentially all within the shaded 
envelope means that this method gives errors that are consistent 
with the analytical errors measured using the data, justifying our 
nominal choice for calculating the errors at low frequencies. 

However, comparing the black lines with the coloured lines 
clearly indicates that in the data and simulations obtaining the 
error bars by measuring the scatter within the leads to 

an underestimation of the errors by approximately 20 %. Given 



the fact that this difference is observed in both the data and the 
simulations, we exclude any instrumental systematic effect as its 
cause and explain it as being due to noise correlations within 
the ^-bins. Such a correlation is expected, since most of the lens 
reconstruction signal in the /'-range of interest to us comes from 
modes in the CMB map within a relatively narrow range at ( ^ 
1500, and so multipoles in the lens reconstruction are correlated. 
We have also checked that the mask induced i-bm coupling is 
negligible, given the bin width we have chosen, and is always 
smaller than 2 %. 

All of these conclusions remain valid up to 353 GHz. 
However, at 545 and 857 GHz, we see by looking at the er- 
rors measured using simulations (solid black for method 4, solid 
coloured for method 5, and dashed coloured for method 2) that 
the errors deduced from the analytical estimates measured from 
the data are substantially higher than those we measure in the 
simulations. This is easily explained through the fact that we are 
omitting any foreground emission in our simulations. The rel- 
ative contribution of Galactic foreground emission is higher at 
low €, which is expected because the Galactic cirrus emission 
has a steep power spectrum. Overall, the amplitude of this con- 
tribution is also consistent with what is seen in Fig. 5. 

Since the scatter within the ^-bins measured in simulations 
(black dashed line) is about 20 % lower than the data-based ana- 
lytical estimates at all frequencies, we use the data-based analyti- 
cal estimates as the basis for our statistical errors. We could alter- 
natively scale the scatter-determined errors by 20 % and obtain 
consistent results. This approach accounts for the foreground 
emission seen at 545 GHz and 857 GHz, but will in practice ne- 
glect the contribution to the errors from the non-Gaussian part 
of the foregrounds. However, we show in Sect. 5.3.5 that this 
contribution is small enough that we can ignore it. 
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The remaining method to discuss is obtained from the cross- 
spectrum of the reconstructed lensing potential in simulated 
maps with the data-measured temperature maps (method 6, 
coloured dashed-dotted line in Fig. A.l). At 545 and 857 GHz 
the CIB signal is dominant over a large i range, and so the er- 
ror obtained from this method is equal to the "signal" terms in 
Eq. 4, which are the orange and green lines in Fig. 5. These 
two lines make up a significant fraction of the total error and 
provide a reasonable approximation to the true error at high-/'. 
However, at low-^ where Galactic emission is important, and at 
100-353 GHz where the CMB and instrumental noise are the 
largest components, the orange and green curves do not accu- 
rately describe the total error. We can see from Fig. A.l that the 
errors obtained using this method are close to the errors mea- 
sured using the other techniques. However, this method will un- 
derestimate the true errors, since the scatter in the CMB and 
noise components is neglected. 

Note that the results presented in Fig. A. 1 are all computed 
using the 40 % Galaxy mask, but we have checked that they hold 
when using the 20 % and 60 % Galaxy masks (which are dis- 
cussed in detail in Sect. 5.2) and that the results show the appro- 
priate / s ky scaling. 
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